rm(list=ls())
## load libraries
library(readxl)
library(kableExtra)
library(ComplexHeatmap)
library(OlinkAnalyze)
library(factoextra)
library(ggpubr)
library(tidyverse)
library(RColorBrewer)
library(limma)
library(ggrepel)
filter = dplyr::filter
select = dplyr::select
palette_treatment <- c(brewer.pal(3, "Set1"))
col_treatment <- list(treatment = c("Healthy" = "white", "Vehicle" = "black", "Isotype" = "#FFD500", "aCD47" = "#F99C1C", "EGFRvIII CAR" = "#008FD5", "EGFRvIII CAR + aCD47" = "#009B49", "EGFRvIII-SGRP CAR" = "#0054A6", "CD19 CAR"= "#EE1C24", "CD19-SGRP CAR"= "#981B1E"))
col_treatment_dataset <- list(treatment = c("Healthy" = "white", "Vehicle" = "black", "Isotype" = "#FFD500", "aCD47" = "#F99C1C", "EGFRvIII CAR" = "#008FD5", "EGFRvIII CAR + aCD47" = "#009B49", "EGFRvIII-SGRP CAR" = "#0054A6", "CD19 CAR"= "#EE1C24", "CD19-SGRP CAR"= "#981B1E"), dataset = c("nov2022" = "lightgrey", "june2023" = "darkgrey", "bridged" = "#3f3f3f"))
my_boxplot_colors <- scale_fill_manual(breaks = c("Healthy", "Vehicle", "Isotype", "aCD47", "EGFRvIII CAR", "EGFRvIII CAR + aCD47", "EGFRvIII-SGRP CAR", "CD19 CAR", "CD19-SGRP CAR"),
values=c("white", "black", "#FFD500", "#F99C1C", "#008FD5", "#009B49", "#0054A6", "#EE1C24", "#981B1E"))
my_pca_colors <- scale_color_manual(breaks = c("Healthy", "Vehicle", "Isotype", "aCD47", "EGFRvIII CAR", "EGFRvIII CAR + aCD47", "EGFRvIII-SGRP CAR", "CD19 CAR", "CD19-SGRP CAR"),
values=c("#dbdad7", "black", "#FFD500", "#F99C1C", "#008FD5", "#009B49", "#0054A6", "#EE1C24", "#981B1E"))
treatment_levels <- c("Healthy", "Vehicle", "Isotype", "aCD47", "CD19 CAR", "CD19-SGRP CAR", "EGFRvIII CAR", "EGFRvIII CAR + aCD47", "EGFRvIII-SGRP CAR")
## Setting up directories
dataDir <- "data"
rdsDir <- "rds"
plotDir <- "plots"
human tumor cells injected in mice lacking Tcells, Bcells, NKcells, but have microglia and macrophages.
Plasma samples from immuno-compromised mice:
my_NPX_data0A <- read_NPX(file.path(dataDir, "Nov2022 Immuno Oncology_Martins_NPX.xlsx"))
my_NPX_data0B <- read_NPX(file.path(dataDir, "Immuno Oncolog_Martins_NPX.xlsx"))
metadata0 <- read_xlsx(file.path(dataDir, "Sample list TM.xlsx"), sheet = 1)%>%
dplyr::rename(SampleID=Name, treatment=Condition)%>%
arrange(SampleID)%>%
mutate(SampleID=gsub(" ","_",paste0("mouse_", SampleID), fixed=TRUE),
treatment = factor(treatment, levels = treatment_levels))
my_NPX_dataA <- my_NPX_data0A%>%
mutate(SampleID=gsub("B_","", gsub(" ","_",paste0("mouse_", SampleID), fixed=TRUE)))%>%
left_join(metadata0, by="SampleID")%>%
mutate(under_LOD=ifelse(NPX<LOD, "yes", "no"),
dataset="nov2022",
treatment = factor(treatment, levels = treatment_levels))
my_NPX_dataB <- my_NPX_data0B%>%
mutate(SampleID=gsub("B_","", gsub(" ","_",paste0("mouse_", SampleID), fixed=TRUE)))%>%
left_join(metadata0, by="SampleID")%>%
mutate(under_LOD=ifelse(NPX<LOD, "yes", "no"),
dataset="june2023",
treatment = factor(treatment, levels = treatment_levels))
my_NPX_data <- bind_rows(my_NPX_dataA, my_NPX_dataB)
table(metadata0$SampleID %in% my_NPX_data$SampleID)
##
## TRUE
## 54
my_NPX_data %>%
mutate(Panel = gsub("Olink ", "", Panel)) %>%
ggplot(aes(x = NPX, fill = dataset)) +
geom_density(alpha = 0.4) +
olink_fill_discrete(coloroption = c("red", "darkblue")) +
set_plot_theme() +
ggtitle("Before bridging normalization: NPX distribution") +
theme(axis.title.x = element_blank(),
axis.title.y = element_blank(),
strip.text = element_text(size = 16),
legend.title = element_blank(),
legend.position = "top")
#### Extract bridging samples
overlapping_samples <- data.frame(SampleID = intersect(my_NPX_dataA$SampleID, my_NPX_dataB$SampleID)) %>%
filter(!str_detect(SampleID, "CONTROL_SAMPLE")) %>% #Remove control samples
pull(SampleID)
npx_before_br <- my_NPX_dataA %>%
dplyr::filter(!str_detect(SampleID, "CONTROL_SAMPLE")) %>% #Remove control samples
dplyr::mutate(Type = if_else(SampleID %in% overlapping_samples,
paste0("nov2022 Bridge"),
paste0("nov2022 Sample"))) %>%
rbind({
my_NPX_dataB %>%
filter(!str_detect(SampleID, "CONTROL_SAMPLE")) %>% #Remove control samples %>%
mutate(Type = if_else(SampleID %in% overlapping_samples,
paste0("june2023 Bridge"),
paste0("june2023 Sample"))) %>%
mutate(SampleID = if_else(SampleID %in% overlapping_samples,
paste0(SampleID, "_new"),
SampleID))
})
md_before_br <- npx_before_br%>%
select(SampleID, treatment, dataset, Type, `Color code`)%>%
unique()%>%
arrange(SampleID)
NPX_noNA <-npx_before_br%>%
dplyr::select(SampleID, Assay, NPX)%>%
spread(Assay,NPX)%>%
select_if(~ !any(is.na(.)))%>%
arrange(SampleID)%>%
column_to_rownames("SampleID")
all.equal(rownames(NPX_noNA), md_before_br$SampleID)
## [1] TRUE
pca <- prcomp(NPX_noNA, scale = T)
fviz_eig(pca)
fviz_pca_ind(pca, axes = c(1, 2), label="none", habillage = md_before_br$Type,
addEllipses=TRUE, ellipse.level=0.95, pointsize = 2,
invisible="quali")
NPX_noNA1 <-NPX_noNA[rownames(NPX_noNA) %in% filter(npx_before_br, dataset=="nov2022")$SampleID,]
md_before_br1 <- md_before_br%>%
filter(dataset=="nov2022")
all.equal(rownames(NPX_noNA1), md_before_br1$SampleID)
## [1] TRUE
pca1 <- prcomp(NPX_noNA1, scale = T) # x is a numeric matrix (row names = patient, colnames= molecule)
fviz_eig(pca1)
fviz_pca_ind(pca1, axes = c(1, 2), label="none", habillage = md_before_br1$Type,
addEllipses=TRUE, ellipse.level=0.95, pointsize = 2,
invisible="quali")
NPX_noNA2 <-NPX_noNA[rownames(NPX_noNA) %in% filter(npx_before_br, dataset=="june2023")$SampleID,]
md_before_br2 <- md_before_br%>%
filter(dataset=="june2023")
all.equal(rownames(NPX_noNA2), md_before_br2$SampleID)
## [1] TRUE
pca2 <- prcomp(NPX_noNA2, scale = T) # x is a numeric matrix (row names = patient, colnames= molecule)
fviz_eig(pca2)
fviz_pca_ind(pca2, axes = c(1, 2), label="none", habillage = md_before_br2$Type,
addEllipses=TRUE, ellipse.level=0.95, pointsize = 2,
invisible="quali")
npx_before_br %>%
mutate(Panel = gsub("Olink ", "", Panel)) %>%
ggplot(aes(x = NPX, group = SampleID, fill=dataset)) +
geom_density(alpha = 0.4) +
set_plot_theme() +
ggtitle("Before bridging normalization: NPX distribution") +
theme(axis.title.x = element_blank(),
axis.title.y = element_blank(),
strip.text = element_text(size = 16),
legend.title = element_blank(),
legend.position = "top")
reference dataset = november2022
overlap_samples_list <- list("DF1" = overlapping_samples,
"DF2" = overlapping_samples)
# Perform Bridging normalization
npx_br_data <- olink_normalization_bridge(project_1_df = my_NPX_dataA,
project_2_df = my_NPX_dataB,
bridge_samples = overlap_samples_list,
project_1_name = "nov2022",
project_2_name = "june2023",
project_ref_name = "nov2022")
# dplyr::glimpse(npx_br_data)
npx_br_data %>%
mutate(Panel = gsub("Olink ", "", Panel)) %>%
ggplot2::ggplot(ggplot2::aes(x = NPX, fill = dataset)) +
ggplot2::geom_density(alpha = 0.4) +
olink_fill_discrete(coloroption = c("red", "darkblue")) +
set_plot_theme() +
ggplot2::ggtitle("After bridging normalization: NPX distribution") +
theme(axis.title.x = element_blank(),
axis.title.y = element_blank(),
strip.text = element_text(size = 16),
legend.title = element_blank(),
legend.position = "top")
npx_after_br <- npx_br_data %>%
dplyr::mutate(Type = ifelse(SampleID %in% overlapping_samples,
paste(Project, "Bridge"),
paste(Project, "Sample")))%>%
dplyr:::mutate(SampleID.orig = SampleID, dataset,
SampleID = paste0(SampleID, "_", dataset))
NPX_noNA <- npx_after_br%>%
dplyr::select(SampleID, Assay, NPX)%>%
spread(Assay, NPX)%>%
select_if(~ !any(is.na(.)))%>%
arrange(SampleID)%>%
column_to_rownames("SampleID")
md_after_br <- npx_after_br%>%
select(SampleID, SampleID.orig, treatment, dataset, Type, `Color code`)%>%
unique()%>%
arrange(SampleID)
NPX <-npx_after_br%>%
dplyr::select(SampleID, Assay, NPX)%>%
spread(Assay,NPX)%>%
arrange(SampleID)%>%
column_to_rownames(var = "SampleID")
all.equal(rownames(NPX), md_after_br$SampleID)
## [1] TRUE
plotDensities(t(NPX), legend=F)
plotDensities(t(NPX), legend = "topright", group= md_after_br$treatment)
i.e: use mean of bridged samples
npx_after_br_mean <- npx_after_br%>%
group_by(SampleID.orig, Assay)%>%
mutate(NPX.ogig=NPX,
NPX=mean(NPX))%>%
ungroup()%>%
distinct(SampleID.orig, Assay, .keep_all = TRUE)%>%
mutate(SampleID=ifelse(SampleID.orig %in% overlapping_samples, SampleID.orig, SampleID),
dataset=ifelse(SampleID.orig %in% overlapping_samples, "bridged", dataset))
NPX_mean <-npx_after_br_mean%>%
dplyr::select(SampleID, Assay, NPX)%>%
spread(Assay,NPX)%>%
arrange(SampleID)%>%
column_to_rownames(var = "SampleID")
md_after_br_mean <- md_after_br%>%
mutate(SampleID=ifelse(SampleID.orig %in% overlapping_samples, SampleID.orig, SampleID),
dataset=ifelse(SampleID.orig %in% overlapping_samples, "bridged", dataset))%>%
distinct(SampleID.orig, .keep_all = TRUE)
plotDensities(t(NPX_mean), legend=F)
plotDensities(t(NPX_mean), legend = "topright", group= md_after_br_mean$treatment)
NPX.norm_mean <- t(normalizeBetweenArrays(t(NPX_mean), method = "cyclicloess"))
all.equal(rownames(NPX.norm_mean), md_after_br_mean$SampleID)
## [1] TRUE
plotDensities(t(NPX.norm_mean), legend=F)
plotDensities(t(NPX.norm_mean), legend = "topright", group= md_after_br_mean$treatment)
NPX.norm_mean.lf <- data.frame(NPX.norm_mean, check.names = F)%>%
rownames_to_column("SampleID")%>%
pivot_longer(values_to = "NPX.norm", names_to = "Assay", !c(SampleID))
npx_after_br_mean <- npx_after_br_mean%>%
left_join(NPX.norm_mean.lf, by=c("SampleID", "Assay"))%>%
relocate(NPX.norm, .after=NPX)
rm(NPX.norm_mean.lf)
# annotation
lbls_col = gsub("_2$","", md_after_br_mean$SampleID.orig)
column_ha = HeatmapAnnotation(
treatment = md_after_br_mean$treatment,
dataset = md_after_br_mean$dataset,
col = col_treatment_dataset, annotation_name_side="left")
## annotation missing freq
miss.per.dataset_mean<-npx_after_br_mean%>%
group_by(dataset, Assay)%>%
summarise(MissingFreq=MissingFreq)%>%
unique()%>%
pivot_wider(names_from = dataset, values_from = MissingFreq)%>%
type.convert()
all.equal(miss.per.dataset_mean$Assay, colnames(NPX.norm_mean))
## [1] "5 string mismatches"
miss.per.dataset_mean<-miss.per.dataset_mean[match(colnames(NPX.norm_mean), miss.per.dataset_mean$Assay),]
all.equal(miss.per.dataset_mean$Assay, colnames(NPX.norm_mean))
## [1] TRUE
row_ha = rowAnnotation(
MissingFreq_june2023 = anno_barplot(miss.per.dataset_mean$june2023, gp = gpar(fill = col_treatment_dataset$dataset["june2023"])),
MissingFreq_nov2022 = anno_barplot(miss.per.dataset_mean$nov2022, gp = gpar(fill = col_treatment_dataset$dataset["nov2022"])), annotation_name_side="top")
# heatmap
hm <- Heatmap(t(scale(NPX.norm_mean, center = T)),
column_labels = lbls_col,
bottom_annotation = column_ha,
right_annotation = row_ha,
na_col = "black",
# col = col_fun,
name = "Z-Score",
column_names_gp = grid::gpar(fontsize = 10),
row_names_gp = grid::gpar(fontsize = 7))
hm
pdf(file.path(plotDir, "Supp_Fig5d_Normalized_IO_prot_in_plasma.pdf"), height =10, width=10)
hm
dev.off()
## quartz_off_screen
## 2
remove NAs
NPX.norm_noNA_mean <- as.data.frame(NPX.norm_mean)%>%
select_if(~ !any(is.na(.)))
all.equal(rownames(NPX.norm_noNA_mean), md_after_br_mean$SampleID)
## [1] TRUE
pca <- prcomp(NPX.norm_noNA_mean, scale = T)
summary(pca)
## Importance of components:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7
## Standard deviation 3.9311 3.1320 2.55332 2.10463 1.96399 1.8687 1.72941
## Proportion of Variance 0.1717 0.1090 0.07244 0.04922 0.04286 0.0388 0.03323
## Cumulative Proportion 0.1717 0.2807 0.35313 0.40235 0.44521 0.4840 0.51724
## PC8 PC9 PC10 PC11 PC12 PC13 PC14
## Standard deviation 1.69400 1.68393 1.62095 1.55213 1.52804 1.46194 1.43134
## Proportion of Variance 0.03188 0.03151 0.02919 0.02677 0.02594 0.02375 0.02276
## Cumulative Proportion 0.54912 0.58063 0.60982 0.63659 0.66254 0.68628 0.70905
## PC15 PC16 PC17 PC18 PC19 PC20 PC21
## Standard deviation 1.36434 1.32921 1.28198 1.26903 1.25160 1.1503 1.14732
## Proportion of Variance 0.02068 0.01963 0.01826 0.01789 0.01741 0.0147 0.01463
## Cumulative Proportion 0.72973 0.74936 0.76762 0.78551 0.80292 0.8176 0.83225
## PC22 PC23 PC24 PC25 PC26 PC27 PC28
## Standard deviation 1.10338 1.05048 1.02724 0.97778 0.93495 0.92772 0.89838
## Proportion of Variance 0.01353 0.01226 0.01172 0.01062 0.00971 0.00956 0.00897
## Cumulative Proportion 0.84578 0.85804 0.86976 0.88038 0.89010 0.89966 0.90863
## PC29 PC30 PC31 PC32 PC33 PC34 PC35
## Standard deviation 0.87035 0.85432 0.81265 0.80314 0.76918 0.71921 0.71486
## Proportion of Variance 0.00842 0.00811 0.00734 0.00717 0.00657 0.00575 0.00568
## Cumulative Proportion 0.91704 0.92515 0.93249 0.93966 0.94623 0.95198 0.95766
## PC36 PC37 PC38 PC39 PC40 PC41 PC42
## Standard deviation 0.66055 0.64295 0.60276 0.5771 0.54667 0.53153 0.52491
## Proportion of Variance 0.00485 0.00459 0.00404 0.0037 0.00332 0.00314 0.00306
## Cumulative Proportion 0.96251 0.96710 0.97114 0.9748 0.97816 0.98130 0.98436
## PC43 PC44 PC45 PC46 PC47 PC48 PC49
## Standard deviation 0.48292 0.45429 0.42070 0.40445 0.37370 0.33992 0.31297
## Proportion of Variance 0.00259 0.00229 0.00197 0.00182 0.00155 0.00128 0.00109
## Cumulative Proportion 0.98695 0.98924 0.99121 0.99303 0.99458 0.99586 0.99695
## PC50 PC51 PC52 PC53 PC54
## Standard deviation 0.30463 0.28030 0.23838 0.21528 7.382e-15
## Proportion of Variance 0.00103 0.00087 0.00063 0.00051 0.000e+00
## Cumulative Proportion 0.99798 0.99885 0.99949 1.00000 1.000e+00
fviz_eig(pca)
fviz_pca_ind(pca, axes = c(1, 2), label="none", habillage = md_after_br_mean$treatment,
addEllipses=TRUE, ellipse.level=0.95, pointsize = 2, invisible="quali")+
my_pca_colors+
my_boxplot_colors
pca1 <- as.data.frame(pca$x)%>%
rownames_to_column("SampleID")%>%
left_join(md_after_br_mean)
pdf(file.path(plotDir, "Fig4a_PCA_plasma_markers_distrib.pdf"), width=10, height = 6)
fviz_pca_ind(pca, axes = c(1, 2), label="none", habillage = md_after_br_mean$treatment,
addEllipses=TRUE, ellipse.level=0.80, pointsize = 2, invisible="quali")+
my_pca_colors+
my_boxplot_colors
ggplot(pca1, aes(PC1,PC2))+
theme_bw()+
geom_point(aes(color=treatment, shape=treatment), size=3)+
stat_ellipse(aes(color=treatment, fill=treatment), geom = "polygon", alpha=0.2)+
geom_hline(yintercept=0, linetype="dashed")+
geom_vline(xintercept=0, linetype="dashed")+
my_pca_colors+
my_boxplot_colors+
scale_shape_manual(breaks = c("Healthy", "Vehicle", "Isotype", "aCD47", "EGFRvIII CAR", "EGFRvIII CAR + aCD47", "EGFRvIII-SGRP CAR", "CD19 CAR", "CD19-SGRP CAR"), values=c(16, 17, 18, 0, 6, 8, 5, 7, 15))
dev.off()
## quartz_off_screen
## 2
all.equal(md_after_br_mean$SampleID, rownames(NPX.norm_mean))
## [1] TRUE
moma <- model.matrix(~ 0 + treatment + dataset, data=md_after_br_mean)
colnames(moma) <- gsub("[^[:alnum:]]+", ".", gsub("treatment", "", colnames(moma), fixed = T))
## Contrast matrix
contrasts.matrix <- makeContrasts(
EGFRvIII.CAR_vs_EGFRvIII.CAR.aCD47 = EGFRvIII.CAR - EGFRvIII.CAR.aCD47,
EGFRvIII.CAR_vs_EGFRvIII.SGRP.CAR = EGFRvIII.CAR - EGFRvIII.SGRP.CAR,
# EGFRvIII.CAR_vs_CD19.SGRP.CAR = EGFRvIII.CAR - CD19.SGRP.CAR,
EGFRvIII.CAR_vs_CD19.CAR = EGFRvIII.CAR - CD19.CAR,
CD19.SGRP.CAR_vs_EGFRvIII.SGRP.CAR = CD19.SGRP.CAR - EGFRvIII.SGRP.CAR,
CD19.CAR_vs_CD19.SGRP.CAR = CD19.CAR - CD19.SGRP.CAR,
EGFRvIII.CAR.aCD47_vs_EGFRvIII.SGRP.CAR = EGFRvIII.CAR.aCD47 - EGFRvIII.SGRP.CAR,
aCD47_vs_EGFRvIII.CAR.aCD47 = aCD47 - EGFRvIII.CAR.aCD47,
levels=moma)
## limma-trend
fit <- lmFit(t(NPX.norm_mean), moma)
fit2 <- contrasts.fit(fit, contrasts.matrix)
fit2 <- eBayes(fit2, trend=T, robust=F) ## robust=TRUE avoids shrinkage of large variances
table(de <- decideTests(fit2, p.value = 0.1))
##
## -1 0 1
## 41 572 31
results <- topTable(fit2, p.value = 1, number = Inf, coef="EGFRvIII.CAR_vs_EGFRvIII.SGRP.CAR")
kable(results) %>%
kable_styling(position = "center", bootstrap_options = c("striped", "hover")) %>%
scroll_box(width = "800px", height = "400px")
| logFC | AveExpr | t | P.Value | adj.P.Val | B | |
|---|---|---|---|---|---|---|
| CCL3 | -0.9623920 | -0.9667431 | -6.6300760 | 0.0000000 | 0.0000028 | 8.7266751 |
| CD27 | 0.5989700 | 0.0434190 | 4.2629293 | 0.0000968 | 0.0044517 | 0.8723171 |
| IL13 | -1.2176303 | -0.7713638 | -3.5113456 | 0.0009978 | 0.0305987 | -1.3492629 |
| FGF2 | -0.3815010 | -1.4941712 | -2.8468567 | 0.0065331 | 0.1242215 | -3.0987850 |
| IL-1 alpha | -0.7644022 | -2.8846819 | -2.8345475 | 0.0067512 | 0.1242215 | -3.1288765 |
| HGF | 0.5530473 | -0.6376821 | 2.5662891 | 0.0135365 | 0.2002119 | -3.7611730 |
| CD5 | 0.5656122 | -0.2037086 | 2.5191764 | 0.0152335 | 0.2002119 | -3.8674044 |
| MUC-16 | 0.3009491 | -1.6923834 | 2.4529876 | 0.0179448 | 0.2063649 | -4.0141162 |
| ARG1 | -0.7029152 | 0.2162739 | -2.3758325 | 0.0216507 | 0.2213182 | -4.1813242 |
| KLRD1 | 0.3425087 | -0.8553006 | 2.2409556 | 0.0298070 | 0.2504705 | -4.4634884 |
| IL8 | 0.8367005 | 1.1624727 | 2.2389350 | 0.0299476 | 0.2504705 | -4.4676154 |
| CD40 | 0.5446794 | -1.1657998 | 2.1036154 | 0.0408073 | 0.3128557 | -4.7370979 |
| IL18 | -0.3684830 | -1.6520022 | -2.0368869 | 0.0473329 | 0.3349713 | -4.8648831 |
| FASLG | -0.2419810 | -0.7954715 | -1.9414224 | 0.0582340 | 0.3652762 | -5.0416837 |
| HO-1 | -0.8897826 | 1.8787318 | -1.8997029 | 0.0636364 | 0.3652762 | -5.1166867 |
| MCP-2 | 0.3254890 | -0.6541967 | 1.8677018 | 0.0680648 | 0.3652762 | -5.1732729 |
| IL12RB1 | 0.3427914 | -1.2190962 | 1.8417367 | 0.0718477 | 0.3652762 | -5.2185780 |
| TRAIL | -0.2379177 | -1.0874565 | -1.8401513 | 0.0720843 | 0.3652762 | -5.2213265 |
| CD8A | 0.3620962 | -1.4420661 | 1.8181448 | 0.0754375 | 0.3652762 | -5.2592668 |
| PTN | 0.4286965 | -2.2701437 | 1.7609400 | 0.0847733 | 0.3899571 | -5.3560339 |
| CASP-8 | -0.5317449 | -1.3645874 | -1.6785476 | 0.0998980 | 0.4373142 | -5.4906354 |
| TWEAK | 0.3084861 | 1.9511467 | 1.6551612 | 0.1045751 | 0.4373142 | -5.5278021 |
| TNF | -0.2953175 | -0.7343338 | -1.5749503 | 0.1219961 | 0.4453495 | -5.6517413 |
| IL15 | 0.6558730 | -0.7330178 | 1.5524183 | 0.1277490 | 0.4453495 | -5.6557620 |
| MMP7 | 0.6150454 | 1.6256308 | 1.5454253 | 0.1289709 | 0.4453495 | -5.6959712 |
| KIR3DL1 | 0.3366558 | -2.5137428 | 1.5405396 | 0.1301554 | 0.4453495 | -5.7032174 |
| CRTAM | -0.2660694 | -1.2714628 | -1.5383033 | 0.1307004 | 0.4453495 | -5.7065273 |
| Gal-9 | 0.1878269 | -1.1532767 | 1.4830529 | 0.1447554 | 0.4756248 | -5.7869183 |
| TNFRSF21 | 0.2100938 | 5.6156598 | 1.4312399 | 0.1589974 | 0.5044055 | -5.8598749 |
| IL10 | 0.2432414 | -1.5527216 | 1.3751041 | 0.1756385 | 0.5336925 | -5.9362337 |
| MMP12 | -0.1643061 | -1.3175578 | -1.3616129 | 0.1798312 | 0.5336925 | -5.9541659 |
| CD4 | 0.2562845 | -1.9008701 | 1.3349838 | 0.1883315 | 0.5414531 | -5.9890811 |
| PDGF subunit B | -0.6601414 | 6.5580033 | -1.3105465 | 0.1963983 | 0.5452567 | -6.0205603 |
| LAP TGF-beta-1 | -0.5032287 | 3.4765154 | -1.2875889 | 0.2042122 | 0.5452567 | -6.0496413 |
| MCP-3 | 0.2168590 | -1.7680756 | 1.2783134 | 0.2074346 | 0.5452567 | -6.0612551 |
| ICOSLG | -0.1447478 | -1.3615353 | -1.2372297 | 0.2221661 | 0.5677579 | -6.1117531 |
| NOS3 | -0.3621967 | -1.5350210 | -1.1761168 | 0.2454888 | 0.6104047 | -6.1840076 |
| CCL23 | -0.1538551 | -0.4200597 | -1.1272682 | 0.2653685 | 0.6424710 | -6.2392806 |
| MIC-A/B | 0.1944967 | -1.3480782 | 1.0539194 | 0.2973240 | 0.6882977 | -6.3180982 |
| DCN | 0.1432047 | -2.1245518 | 1.0324737 | 0.3071498 | 0.6882977 | -6.3401889 |
| CXCL13 | -0.2062826 | -2.0376865 | -1.0183412 | 0.3137450 | 0.6882977 | -6.3545094 |
| VEGFA | 0.1707794 | -1.4025570 | 1.0173249 | 0.3142229 | 0.6882977 | -6.3555319 |
| IL33 | 0.1264109 | -2.2348357 | 0.9458131 | 0.3490957 | 0.7302897 | -6.4250240 |
| CXCL10 | -0.2895126 | -2.1843327 | -0.9318619 | 0.3561845 | 0.7302897 | -6.4380139 |
| IFN-gamma | -0.7195532 | 4.2931944 | -0.9193190 | 0.3627392 | 0.7302897 | -6.4419794 |
| CXCL12 | 0.1666708 | -0.3822878 | 0.9144832 | 0.3651448 | 0.7302897 | -6.4539349 |
| ANGPT2 | -0.1222052 | -1.9008921 | -0.8919054 | 0.3770009 | 0.7379592 | -6.4741869 |
| GZMH | -0.1335740 | -1.6666856 | -0.7576131 | 0.4524758 | 0.8530743 | -6.5844914 |
| CCL19 | -0.1248333 | -2.0365889 | -0.7510865 | 0.4563561 | 0.8530743 | -6.5894066 |
| GZMA | 0.1848707 | -2.2433253 | 0.7389426 | 0.4636273 | 0.8530743 | -6.5984417 |
| GZMB | -0.0884843 | -1.4916136 | -0.6530481 | 0.5169167 | 0.8790867 | -6.6582309 |
| CCL17 | 0.2417116 | 1.1842148 | 0.6527681 | 0.5170957 | 0.8790867 | -6.6584140 |
| TNFSF14 | 0.1735614 | 0.6228065 | 0.6264163 | 0.5340818 | 0.8790867 | -6.6752981 |
| IL7 | 0.1096035 | -0.0869920 | 0.6234895 | 0.5359862 | 0.8790867 | -6.6771311 |
| PDCD1 | 0.1895465 | -0.8762071 | 0.5924066 | 0.5564263 | 0.8790867 | -6.6960772 |
| TIE2 | 0.1008270 | -1.8918147 | 0.5780468 | 0.5660002 | 0.8790867 | -6.7045077 |
| CCL4 | -0.0872828 | -0.5336423 | -0.5757577 | 0.5675339 | 0.8790867 | -6.7058328 |
| IL12 | 0.1883311 | -1.6968559 | 0.5579164 | 0.5795574 | 0.8790867 | -6.7159828 |
| IL5 | -0.1798104 | -1.4721251 | -0.5523870 | 0.5833086 | 0.8790867 | -6.7190646 |
| CSF-1 | -0.1548806 | -1.3511787 | -0.5491891 | 0.5854834 | 0.8790867 | -6.7208330 |
| CAIX | -0.0781767 | -1.1660219 | -0.5477622 | 0.5864551 | 0.8790867 | -6.7216189 |
| CD40-L | 0.1030995 | -2.1076827 | 0.4995181 | 0.6197518 | 0.8790867 | -6.7469989 |
| TNFRSF9 | -0.1107819 | -0.2974723 | -0.4969672 | 0.6215359 | 0.8790867 | -6.7482765 |
| Gal-1 | -0.1074847 | 0.2941486 | -0.4936574 | 0.6238541 | 0.8790867 | -6.7499246 |
| CXCL1 | 0.1070344 | -0.6135119 | 0.4870278 | 0.6285091 | 0.8790867 | -6.7531928 |
| PGF | 0.0868608 | -0.6154360 | 0.4785671 | 0.6344721 | 0.8790867 | -6.7573003 |
| TNFRSF4 | -0.1413895 | -1.9483573 | -0.4530296 | 0.6526182 | 0.8790867 | -6.7692657 |
| CD83 | 0.0599345 | -0.7571030 | 0.4504176 | 0.6544865 | 0.8790867 | -6.7704529 |
| CXCL5 | -0.0559320 | -2.1444526 | -0.4436813 | 0.6593150 | 0.8790867 | -6.7734833 |
| ANGPT1 | 0.0500481 | 8.2381087 | 0.4180397 | 0.6778277 | 0.8908592 | -6.7846038 |
| LAMP3 | 0.0666158 | -1.2008492 | 0.4047024 | 0.6875378 | 0.8908940 | -6.7901284 |
| CD28 | 0.0575505 | -1.0302463 | 0.3671227 | 0.7151803 | 0.9138414 | -6.8047376 |
| CCL20 | 0.0732618 | -1.1730268 | 0.3229705 | 0.7481544 | 0.9270812 | -6.8200942 |
| IL6 | 0.1885328 | -0.3982624 | 0.3207243 | 0.7498453 | 0.9270812 | -6.8208231 |
| EGF | 0.1112740 | -0.1910594 | 0.3128634 | 0.7557727 | 0.9270812 | -6.8233344 |
| ADA | 0.0582176 | -2.2622091 | 0.2921652 | 0.7714500 | 0.9338605 | -6.8296498 |
| IL4 | 0.0915919 | -1.1894395 | 0.2591126 | 0.7966833 | 0.9395703 | -6.8388419 |
| LAG3 | 0.0509785 | -2.3299600 | 0.2392899 | 0.8119243 | 0.9395703 | -6.8438274 |
| NCR1 | 0.0453799 | -1.6939142 | 0.2274449 | 0.8210673 | 0.9395703 | -6.8466177 |
| CXCL11 | 0.1011961 | -0.4275654 | 0.2181226 | 0.8282807 | 0.9395703 | -6.8487143 |
| MCP-1 | 0.0870971 | 0.7340598 | 0.2158379 | 0.8300509 | 0.9395703 | -6.8492148 |
| MCP-4 | -0.0378219 | 0.9295177 | -0.2022974 | 0.8405600 | 0.9395703 | -6.8520730 |
| IL2 | -0.0679478 | -1.8765798 | -0.1917905 | 0.8487350 | 0.9395703 | -6.8541636 |
| CD70 | -0.0217853 | -0.3738565 | -0.1800771 | 0.8578686 | 0.9395703 | -6.8563629 |
| PD-L1 | -0.0669698 | -1.3769988 | -0.1578187 | 0.8752778 | 0.9449407 | -6.8601610 |
| VEGFR-2 | -0.0216313 | -1.4390177 | -0.1423947 | 0.8873788 | 0.9449407 | -6.8624996 |
| CD244 | -0.0191346 | -1.9857875 | -0.1290871 | 0.8978412 | 0.9449407 | -6.8643244 |
| PD-L2 | 0.0112252 | -1.8406661 | 0.1153445 | 0.9086649 | 0.9449407 | -6.8660213 |
| CX3CL1 | -0.0138398 | -1.3156902 | -0.0984455 | 0.9219985 | 0.9449407 | -6.8678464 |
| CXCL9 | 0.0107076 | -1.3096721 | 0.0954068 | 0.9243986 | 0.9449407 | -6.8681440 |
| ADGRG1 | 0.0117546 | -1.5616801 | 0.0696628 | 0.9447585 | 0.9551405 | -6.8702910 |
| TNFRSF12A | 0.0098702 | 4.7079791 | 0.0386012 | 0.9693723 | 0.9693723 | -6.8719901 |
# Categorize Results based on P-value & FDR for plotting
results$Color <- "NS or FC < 0.5"
results$Color[results$P.Value < 0.05] <- "P < 0.05"
results$Color[results$adj.P.Val < 0.05] <- "adj.P.Val < 0.05"
results$Color[results$adj.P.Val < 0.001] <- "adj.P.Val < 0.001"
results$Color[abs(results$logFC) < 0.5] <- "NS or FC < 0.5"
results$Color <- factor(results$Color,
levels = c("NS or FC < 0.5", "P < 0.05",
"adj.P.Val < 0.05", "adj.P.Val < 0.001"))
p1 <- ggplot(results,
aes(x = logFC, y = -log10(P.Value),
color = Color, label = rownames(results))) +
geom_vline(xintercept = c(0.5, -0.5), lty = "dashed") +
geom_hline(yintercept = -log10(0.05), lty = "dashed") +
geom_point() +
labs(x = "Enriched in EGFRvIII-SGRP <- log2(FC) -> Enriched in EGFRvIII",
y = "Significance, -log10(P)",
color = "Significance") +
scale_color_manual(values = c(`adj.P.Val < 0.001` = "dodgerblue",
`adj.P.Val < 0.05` = "lightblue",
`P < 0.05` = "orange2",
`NS or FC < 0.5` = "gray"),
guide = guide_legend(override.aes = list(size = 4))) +
scale_y_continuous(expand = expansion(mult = c(0,0.05))) +
geom_text_repel(min.segment.length = 0.3) +
theme_bw(base_size = 16) +
theme(legend.position = "bottom")
p1
pdf(file.path(plotDir, "Fig4b_volcano_EGFRvIII_CAR_vs_EGFRvIII_SGRP_CAR.pdf"))
p1
dev.off()
## quartz_off_screen
## 2
interest.prot2 <- unique(rownames(topTable(fit2, p.value = 0.05, number = Inf, coef=c(1:7))))
comparisons <- list(c("EGFRvIII CAR", "EGFRvIII CAR + aCD47"),
c("EGFRvIII CAR", "EGFRvIII-SGRP CAR"),
c("EGFRvIII CAR", "CD19 CAR"),
c("CD19-SGRP CAR", "EGFRvIII-SGRP CAR"),
c("CD19 CAR", "CD19-SGRP CAR"),
c("EGFRvIII CAR + aCD47", "EGFRvIII-SGRP CAR"),
c("aCD47", "EGFRvIII CAR + aCD47"))
data <- npx_after_br_mean%>%
filter(Assay %in% interest.prot2,
!is.na(NPX))%>%
mutate(Assay=factor(Assay, levels = interest.prot2),
treatment=factor(treatment, levels = treatment_levels))
stat.test <- data%>%
group_by(Assay)%>%
rstatix::wilcox_test(NPX.norm~treatment,
comparisons = comparisons)%>%
rstatix::add_xy_position()%>%
mutate(grp=paste0(group1, "_", group2))%>%
group_by(grp)%>%
mutate(p.adj=p.adjust(p, method = "BH", n = length(p)),
p.adj.signif=ifelse(p.adj<0.001, "***", ifelse(p.adj<0.01, "**", ifelse(p.adj<=0.05, "*", "ns"))))%>%
ungroup()
write.csv(select(stat.test, -groups), file.path(plotDir,"stat_test.csv"), row.names = F )
p2 <- ggplot(data)+
theme_bw()+
geom_boxplot(aes(x=treatment, y=NPX.norm, fill=treatment), outlier.shape = NA, alpha=0.7)+
my_boxplot_colors+
geom_jitter(aes(x=treatment, y=NPX.norm, color=under_LOD), alpha=0.8)+
scale_color_manual(values=c("black", "#999999"))+
theme(axis.text.x=element_blank(), #remove x axis labels
axis.ticks.x=element_blank())+ #remove x axis ticks
ggtitle("normalised values")+
facet_wrap(~Assay, scales="free", ncol=3)
p2
pdf(file.path(plotDir, "Fig4c_norm_exp_DE_plasma_markers.pdf"), width = 11, height = 13)
p2
dev.off()
## quartz_off_screen
## 2
p3 <- ggplot(filter(npx_after_br_mean, Assay=="IL6"))+
theme_bw()+
geom_boxplot(aes(x=treatment, y=NPX.norm, fill=treatment), outlier.shape = NA, alpha=0.7)+
my_boxplot_colors+
geom_jitter(aes(x=treatment, y=NPX.norm, color=under_LOD), alpha=0.8)+
scale_color_manual(values=c("black", "#999999"))+
theme(axis.text.x=element_blank(), #remove x axis labels
axis.ticks.x=element_blank())+ #remove x axis ticks
ggtitle("IL6 normalised values")
p3
pdf(file.path(plotDir, "Supp_Fig7c_Plasma_IL6.pdf"))
p3
dev.off()
## quartz_off_screen
## 2